function haversine_f(deglat1,deglon1,deglat2,deglon2) result (dist)

    use iasing_osse
   
    implicit none
    ! Great circle distance 
    real(kind=JPR4),intent(in) :: deglat1,deglon1,deglat2,deglon2
    real(kind=JPR4)            :: a,c,dist,dlat,dlon,lat1,lat2
   
    dlat = real(DEG2RAD,4) * (deglat2-deglat1)
    dlon = real(DEG2RAD,4) * (deglon2-deglon1)
    lat1 = real(DEG2RAD,4) * (deglat1)
    lat2 = real(DEG2RAD,4) * (deglat2)
    a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2
    c = 2*atan2(sqrt(a),sqrt(1-a))
    dist = real(RADEARTH,4) * c
    
end function haversine_f
